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Abstract. 

The equation governing the streaming of a quantity down its gradient superfieially looks similar 
to the simple constant velocity advection equation. In fact, it is the same as an advection equation 
if there are no local extrema in the computational domain or at the boundary. However, in general 
when there are local extrema in the computational domain it is a non-trivial nonlinear equation. 
The standard upwind time evolution with a CFL-limited time step results in spurious oscillations 
at the grid scale. These oscillations, which originate at the extrema, propagate throughout the 
computational domain and are undamped even at late times. These oscillations arise because of 
unphysically large fluxes leaving (entering) the maxima (minima) with the standard CFL-limited 
explicit methods. Regularization of the equation shows that it is diffusive at the extrema; because 
of this, an explicit method for the regularized equation with A.t (x A.x^ behaves fine. We show that 
the implicit methods show stable and converging results with At oc Ax; however, surprisingly, even 
implicit methods are not stable with large enough timesteps. In addition to these subtleties in the 
numerical implementation, the solutions to the streaming equation are quite novel: non-differentiable 
solutions emerge from initially smooth profiles; the solutions show transport over large length scales, 
e.g., in form of tails. The fluid model for cosmic rays interacting with a thermal plasma (valid at 
space scales much larger than the cosmic ray Larmor radius) is similar to the equation for streaming 
of a quantity down its gradient, so our method will find applications in fluid modeling of cosmic rays. 
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1. Introduction. Cosmic rays (energetic particles moving close to the speed of 
light) are an important dynamic and energetic component of our Galaxy. The average 
cosmic ray pressure in our Galaxy is comparable to the kinetic, thermal, and magnetic 
pressures [2] . It is important to study the interaction of cosmic rays with the thermal 
plasma to understand the dynamics and the energy balance in our Galaxy and in 
other astrophysical objects, e.g., clusters of galaxies. 

Cosmic rays are collisionless, with the Coulomb collision time ^ the dynamical 
timescales. However, cosmic rays do not just leave the system at the speed of light. 
Cosmic rays are effectively coupled to the thermal plasma by self-generated Alfven 
waves that scatter them (e.g., [1]; efficient scattering occurs only when the cosmic 
rays have a sufficiently large number density; this is usually the case for the majority 
of cosmic ray particles in the Galaxy; e.g., [M]). These waves are generated by a fast, 
self-generated instability that arises when the cosmic ray streaming velocity (velocity 
of the cosmic ray fluid relative to the thermal plasma) exceeds the local Alfven speed. 
As a result of this instability, the cosmic ray fluid streams at the local Alfven speed 
(down the cosmic ray pressure gradient) with respect to the thermal plasma. Other 
mechanisms, like MHD turbulence (e.g., [15]) and non-resonant instabilities (e.g., [1]) 
can also scatter cosmic rays but the qualitative picture of streaming relative to the 
thermal plasma should hold even in these cases. For simplicity, we focus on the self- 
generated streaming instability and assume that the cosmic ray fluid streams relative 
to the thermal plasma at the local Alfven speed. 
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Although a fuUy kinetic description is necessary to study the excitation of Alfven 
waves (at the cosmic ray Larnior radius scale) and their damping on the thermal 
plasma, a fluid description, where the cosmic ray fluid streams at the local Alfven 
speed relative to the thermal plasma, suffices at scales much larger than the cosmic 
ray Larmor radius (e.g., the typical cosmic ray Larmor radius is <C the thickness of 
our Galaxy, and hence a fluid model of cosmic rays is sufficient to study the large 
scale dynamics and energetics). The equation governing the cosmic ray pressure is 
(generalization of 1-D equations from 8, ) 

(1-1) -Q^ +^ ■ {Pcu) +V ■ i-pcVs) = -Pc—^ , 

where Pc is the cosmic ray pressure and u is the thermal plasma velocity. The stream- 
ing velocity Vg = — sgn(i?- Vpc)-B/-\/47r/9, where B is the magnetic field vector and p is 
the plasma density, is equal to the Alfven velocity, but down the cosmic ray pressure 
gradient. Equation (jl.ip can be derived by replacing u hy u + Vs in the standard 
internal energy equation for a 7 = 4/3 (relativistic) fluid. The last term on the right 
hand side of Eq. (jl.ip is a loss term since cosmic rays work in pushing the thermal 
plasma as they stream down their own pressure gradient; there is a compensating 
heating term in the plasma internal energy equation. This irreversible heating can 
be understood in the kinetic framework: cosmic ray fluid streaming at speeds slightly 
larger than the local Alfven speed results in (usually small amplitude) Alfven waves 
at the cosmic ray Larmor radius scale. These waves efficiently scatter the cosmic 
ray particles. In addition to elastic scattering, cosmic rays also lose their energy to 
Alfven wave packets because the self-generated waves are moving, on average (in the 
fluid sense), away from the cosmic ray particles (this process is the inverse of Fermi 
acceleration). Thus, a part of cosmic ray energy is converted into Alfven waves which, 
in a steady state, damp on the thermal plasma. In fluid modeling (valid at large time 
and space scales), it is assumed that the cosmic rays directly lose their energy to the 
thermal plasma via the last term in Eq. (|l.ip . 

Eq. (jl.ip can be solved by operator splitting, where each term is updated individ- 
ually. The terms involving u are the usual hyperbolic terms and can be treated with 
the standard methods. The last term on the right hand side can be implemented such 
that it is an exact exchange term between the plasma and cosmic ray internal energy 
equations; this term will probably need to be upwinded. The streaming term (third 
term on the left hand side) is the main subject of this paper, and will be examined 
using the following conservative one-dimensional model equation describing streaming 
down the gradient: 

where f{x,t) is a function of space and time, v — —sgii(df/dx), and sgn(s) — 1 
if s > 0, if s = 0, and — 1 if s < 0. The streaming equation is very similar in 
appearance to the simple advection equation (but very different in nature as we shall 
see), but differs from a variable-coefficient advection equation as the advection speed 
(v) depends, not on x, but on the gradient of /. From here on we focus on the simpler 
Eq. ()1.2p and effective numerical methods to solve it. Once numerical methods for 
solving Eq. (11.21) are available, the streaming part of Eq. (|l.ip can be implemented 
in an analogous fashion. 
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Integrating Eq. (|1.2p in space (from -Ax/2 to Ax/2) and time (from —At/2 to 
At/2) one obtains, 



(1.3) Ax [/o,o+ - /o,o-] + At -/sgn(5//ax)|o+ + fsgn{df/dx)\,^ ^ 



0. 



Evaluating Eq. (jl.3l) at a maximum of / (similar argument applies at a minimum), 
one obtains 

Aa^ [/o,o+ - /o.o-] + At[/o+,o + /o-,o] = 0, 

because the sgn function changes sign across the extremum. Assuming that / is 
continuous in x, this reduces to 

Aa; [/o,o+-/o.o-] +2At/o,o = 0. 

Assuming that / is also continuous in time, this is equivalent to Axdf /dt + 2f = 0. In 
the limit Ax — )■ 0, this implies / = 0, an apparent contradiction. This contradiction 
also applies to the cosmic ray pressure equation (Eq. 11.11 the last term on the right 
hand side of Eq. 11.11 reduces to zero on integrating in space at an extremum) . A 
resolution of this apparent contradiction is that, rather than / = 0, df /dt — )■ —2f/Ax; 
i.e., / changes infinitely fast at an extremum. A similar behavior is observed for the 
diffusion equation when applied at a point where / is continuous but the gradient of / 
changes discontinuously. The analogy with the diffusion equation is deeper; in section 
[3] we show that the regularized equation is diffusive (and hence parabolic) in nature 
at extrema. However, there is a difference. With the diffusion equation, the initial 
discontinuity in df/dx is removed instantly, and df/dt is well behaved at later times. 
For the streaming equation (Eq. 11.21) . however, df/dt —2f/Ax at the maximum 
for all times! This is not mathematically well behaved in the limit Ax — )■ 0. This is 
a problem for both the model equation (Eq. II. 2[) and the cosmic ray equation (Eq. 
II. 1[) . However, when we regularize the sgn(/') function by tanh(/'/e) (/' = df/dx; 
see section[3]for the details of regularization that we use), the flux vf = —f tanh(/'/e) 
instantly becomes continuous (see Fig. 14.51 for the plots of vf at different times) at 
the extrema and df/dt is well behaved, just as it is for the diffusion equation. Thus, 
only in the regularized form can we make sense of the streaming equation, and it is 
reassuring that a unique solution is obtained for the regularized equation (Eq. 13. ip 
in the limit Ax — ?> and e — > 0, for a sufficiently small timestep At . 

The methods for treating the standard hydrodynamic and magnetohydrodynamic 
(MHD) equations are quite mature by now, but entirely unexpected numerical prob- 
lems arise when new physics (as in Eq. II. ip is added. Thus, the numerical implemen- 
tation of these new terms must be tested rigorously. For example, simple centered 
finite differencing of anisotropic thermal conduction can give rise to negative temper- 
atures for large temperature gradients, an unphysical result ([H]). Similarly, here we 
show that a naive implementation of Eq. (|1.2p results in unphysical oscillations. 

The paper is organized as follows. Section ([2]) shows that a standard explicit 
implementation of Eq. (|1.2p results in spurious oscillations of large amplitude. It 
is shown that a much smaller timestep with At cx Ax^ can remove oscillations and 
results in a physically consistent evolution. In section (|3]) we regularize Eq. (|1.2p and 
develop implicit methods based on it. With test problems (with both continuous and 
discontinuous initial profiles) we show that physically consistent solutions are obtained 
using implicit methods with a sufficiently small timestep At oc Ax. In section ([Sj we 
discuss the nature of the solutions, and conclude with the implications of our methods 
for the fiuid modeling of cosmic rays. 
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Fig. 2.1. Profiles at t = (solid black line) and t = 0.02 for the sine wave initial condition 
(Eq. 1 2. j I) with different CFL numbers. Number of grid points is 128. Periodic boundary conditions 
are used. The profiles with cn=0.001 and cn=10~'^ almost overlap. 

2. Non-monotonicity with explicit schemes using the standard CFL 
time-step. The standard methods for solvmg hyperboUc equations fail miserably for 
solving Eq. (|1.2p . as demonstrated by a periodic test problem with a smooth initial 
profile given by 

(2.1) /(x,0) = 1 + 0.5 sin(27rx). 

We show that a standard explicit update scheme with a reasonable Courant-Friedrichs- 
Lewy (CFL) number (Ai = cnAa;/|w|, where cn is the CFL number, \v\ = 1 for Eq. 
I1.2[) gives rise to spurious oscillations even for a smooth initial condition like Eq. (|2.ip . 
Figure (12.11) shows the initial profile and profiles at i = 0.02 for runs with different 
CFL numbers; the number of grid points is 128. Eq. (|1.2p is evolved using an up- 
wind method with a van Leer slope limiter (e.g., see [5]); such oscillations occur for 
all explicit methods. As mentioned in the introduction, the streaming equation (Eq. 
I1.2p is mathematically well-behaved only when the discontinuous sgn(/') function is 
replaced by its smoothened form tanh(/'/e); thus we use this smoothened form with 
e = 0.1 to calculate the fluxes. If we use the sgn function instead of the smoothened 
form, grid scale oscillations in / arise even for an arbitrarily small At (oscillation am- 
plitude becomes smaller with a smaller At though). The profile with cn=0.5 shows 
large amplitude oscillations originating from the initial extrema. Although the ampli- 
tude of oscillations with cn=0.1 is smaller, the oscillations are more spread out; this 
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Fig. 2.2. Zoom in on the profiles at first two timesteps for cn=0.5 in Figure i2. It : the initial 
profile is marked by a solid line and triangles; the profile at the first timestep by dotted line and 
sguares; and the profile at the second timestep by dashed line and circles. Grid-scale oscillations 
originate at the local extrema and propagate throughout the box with time. 

is because the number of timesteps to reach t = 0.02 is larger with a smaher CFL 
number (hence spreading out the osciUations) . As the CFL number is reduced further 
the profiles converge to a reasonable-looking /. Profiles with cn^O.OOl show small 
amplitude oscillations, but the profile with cn=10^^ does not. 

These oscillations originate from the local extrema. With a large timestep (even 
a CFL-limited timestep is large for the extrema), the fluxes (on both sides) leaving 
the initial maximum are big enough to reduce the initial maximum way below its 
initial value (see Fig. 12. 2p . thus giving rise to an unphysical dip at the location of 
the initial maximum. Since the dip is a local minimum now, there are large fluxes 
entering those grid points in the next timestep, again making them a local maximum 
and simultaneously making the adjoining grid cells a minimum. Now two new local 
minima are created in cells adjoining the initial maximum, in addition to the initial 
maximum. In this way oscillations propagate away from the initial extremum as the 
solution is advanced, resulting in oscillations which spread over the whole domain 
with time. 

One can obtain a limit on the timestep such that new local extrema are not 
created. Consider a maximum (at /o) shown schematically in Figure (|2.3p : we consider 
the extrema because the oscillations originate there. Considering a piecewise constant 
reconstruction, let us apply an upwind update of Eq. (jl.2p for a single timestep, 
such that, /o"+^ = - 2At/Axf^, = A" + At/Axif^ ~ f^). Now imposing 

monotonicity after a single timestep, i.e., ff^^^ > fi~^^ gives, 

(2.2) At < |/"|AxV(4|«|/), 

since - w \ f"\Ax^/2 and f^' > \ f"\Ax^, where /" (= d'^f/dx^) is evaluated at 
the extremum. Notice that the timestep restriction will be worse if even higher order 
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Fig. 2.3. A schematic representation of a maximum at /q. Unphysically large fluxes leaving /q 
with the standard CFL-limited timestep result in oscillations seen in Figs. S2. It and \2.S!jl . 

derivatives vanish, and also that |/"| decreases as the initial profile smoothens with 
time. A similar condition for monotonicity is obtained at a minimum. Although we 
have considered an upwind method above, a similar timestep constraint is obtained 
for any explicit method. A result similar to Eq. (|2.2I) is obtained by analyzing the 
regularized form of Eq. (|1.2p in section (jS]). The value of |/"| at extrema for the 
initial profile in Eq. (|2.ip is 27r^ and the limit on timestep given by Eq. (12. 2p is 
At < {tt'^Ax^/3)Ax, which on using Ax = 1/128 gives Ai < 2 x IQ-^Ax. This is 
consistent with oscillations seen in Figure (|2.ip for cn> 0.001 and not for cn^lO^^. 

3. Regularization. Since the advection velocity v can change discontinuously 
at extrema, the differential form in Eq. (|1.2|) is well posed only in its regularized 
form. A similar situation occurs for the Euler/Burger's equations in the absence of 
viscosity. Motivated by the use of explicit viscosity to make Euler/Burger's equa- 
tions well-posed, we introduce a physically motivated regularization of Eq. (|1.2p . A 
straightforward regularization is obtained by replacing sgn(/') in Eq. (|1.2p by its 
smooth approximation tanh(/'/e), where e is a small parameter. 

Thus the regularized form of Eq. (|1.2p becomes, 



side of Eq. p.2p is an advective term, the term on the right hand side of Eq. p.2p is 
a diffusive term, which at an extrcmum (where /' — 0) is ff"/e, a diffusive term with 
the diffusion coefficient //e. There is a physical justification for such a regularization. 
As mentioned before, away from the cosmic ray pressure extrema, cosmic ray fiuid 
streams down its pressure gradient along magnetic field lines; in addition to this 
streaming there is also diffusion (due to pitch angle scattering by the self-excited 
Alfven waves, which is typically small, but should be added on the right hand side of 
Eq. |1.1] ) along magnetic field lines. At cosmic ray pressure extrema, the cosmic ray 
particles have an equal probability to move in either direction (along magnetic field 




Numerical Implementation of Streaming Down the Gradient: Application to Fluid Modeling of Cosmic Rayi 



lines) at the speed of light. However, they are scattered after traveling a distance of 
the order of the Larmor radius, thus resulting in a diffusive behavior at extrema. Just 
as viscosity can be ignored at all locations except close to shocks for Euler/Burger's 
equations, the diffusive behavior of Eq. (13. 2[) is only effective at extrema. 

The time step limit for the stability of the diffusive term in Eq. (|3.2I) treated 
explicitly is 

(3.3) At < Ax^e/2f; 

on taking e ^ ] Axjl? (L is the box-size), this scales in the same way as Eq. (|2.2p . 
Although this restriction on timestep follows from the linear stability analysis of the 
discretized diffusion equation, explicit methods with larger timesteps do not blow 
up because the large diffusion is only concentrated in an infinitesimally small region 
where /" « 0. Away from extrema, the term on the left hand side is like advection. 
Eq. p.2p clearly shows that Eq. (II. 2p is not hyperbolic, and because of that the 
methods for hyperbolic equations do not work for this. Since the regularized form 
in Eq. (j3.ip is more transparent and has a less restrictive timestep requirement, we 
will work with it for testing different explicit and implicit methods. We only focus on 
conservative schemes because the underlying equation (Eq. II. 2|) is conservative. 

3.1. Explicit methods. Explicit methods are not as attractive as implicit meth- 
ods for this problem because of the restrictive timestep (Eq. |3.3| ) for stability, but 
we discuss them for completeness. With the standard CFL timestep, both the Lax- 
Wendroff and upwind methods give very similar large amplitude oscillations, even 
with a smooth initial profile (see Fig. |2.1| ). With a smaller timestep scaling like Acc^ 
(for a fixed e) both methods give physically consistent results. 

3.1.1. Lax-Wendroff method. The second order accurate Lax-Wendroff dis- 
cretization {^) for Eq. p.ip is obtained by observing that 

9/ 



/(x,i-hAt) -/(x,t) + At 



dt 2 9t2 



where partial derivatives (/ = df /dt) are evaluated at (a;, t). Using, f — [f tanh(/'/e)]' 
(Eq. |3.1| ). / = [9(/ tanh(/'/e))/9<]', and expanding the time derivatives in terms of 
space derivatives, one gets, 

^^-^^ = ^ |/tanh(/7e) + — ^tanh(/7e)[/ tanh(/7e)]' + 4ech2(/7e)[/tanh(/7e)]" 
(3.4) 

where quantities on the right hand side are simply averaged to be centered at {i,n). 
The above discretization is clearly conservative. 

3.1.2. Upvirind method. The explicit upwind discretization of Eq. p.ip is 

given by 

where the flux (F^^^y^^) is calculated in an upwind fashion using a standard first order 
slope limiter (we are using the van Leer limiter; |13l The flirx equals 



(/r + a.l^^), for.>0. 



for V < 0, 
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where v — — tanh(/'/e)"_^2/2' '^i ~ (-^ ^ slope limiter like van 

Leer's). 

3.2. Implicit methods. Implicit methods are attractive because longer timesteps, 
scaling as Aa; and not as Ax^, can be used to give stable and converging results as 
the resolution is increased. 

3.2.1. Linearized implicit method. Assuming smoothness of /', we can write 
the backward-Euler expansion of tanh(/'/e) evaluated at timestep (n + 1) in Eq. p.ip 

as 

(3.7) tanh(/7e)"+i = tanh(/7e)" + lsech^(/7e)" - /''") . 

In writing above we have linearized the implicit approximation in The lin- 

earized (semi-)implicit flux for the discretized equation (centered in space), using 
above, can be written as 



(3-8) F-^^ ^ /r+i 



n+1 



tanh(/7e)r+i/2 + -^sech^(/7e)r+i/2/r 

tanh(/7e)iVi/2 - ^sech2(/7e)r+i/2/;;i ) 



which is linear in /"^^. Thus the implicit form of Eq. (|3.1[) . using above, can be 
evolved using the standard tridiagonal solver (since we are imposing periodic boundary 
conditions the tridiagonal solver is combined with the Sherman-Morrison formula; e.g., 

my 

3.2.2. Nonlinear implicit method. Instead of linearizing the implicit equa- 
tion, as we did in the previous section, we can use a fully-implicit nonlinear method. 
The nonlinear implicit equation can be solved using a Krylov subspace method like 
the generalized minimal residual method (GMRES; [10]). Eq. (|3.ip can be discretized 
as a nonlinear matrix equation M(/""'"^)/"+^ = /", where the matrix-vector product 
to be used in GMRES is given by 

(3.9) MV, = V.-^ {V^+l/2 tanh(y7e),;+i/2 - ^,-1/2 tanh(y7e),_i/2} • 

The matrix equation is solved iteratively until the relative error is small enough (we 
choose 10""'^^ as the relative error tolerance). Since matrix M is diagonally dominant, 
for a sufRciently small At, this implicit scheme converges rapidly. 

3.2.3. Semi-implicit method. The implicit methods we consider are formally 
only first order accurate in time, but can be made second order accurate by combining 
implicit and explicit methods, similar to the Crank-Nicolson scheme for the diffusion 
equation. The regularized equation (Eq. f3.1| ) can be discretized as 

Ai " ^ LA+i/2tanh(/7e),+i/2 - /,-i/2 tanh(/7e),_i/2] " 

(3.10) + [/,+i/2 tanh(/7e),+i/2 - /,_i/2 tanh(/7e),_i/2]"+' , 

where the update is split into an explicit and implicit part. We update the explicit 
piece using the Lax-Wendroff method and the implicit piece using the nonlinear im- 
plicit scheme. Results for different test problems are similar as the implicit schemes. 
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(a) At = eAx2/3 



(b) At = 0m2Ax 



Fig. 4.1. The initial profile (Eq. [2771 ) and profiles at t = 0.02 using (a) At = eAx^/3 (e = 0.1) 
and (b) At = 0.002 Ax for explicit schemes. Results with upwind and Lax-Wendroff methods look 
similar. Unlike with At oc Ax^ , with At oc Ax the results do not converge as the spatial resolution 
is increased. 



Somewhat surprisingly, the convergence with increasing number of grid points for the 
semi-implicit {g = 1/2 in Eq. |3.10] ) scheme is not very different from the fully implicit 
methods (see Fig. |4.4j V Notice that /' is not continuous (discontinuity occurs where 
the flattened maxima/minima is connected to the advecting part of the solution; e.g., 
Figs. |4.1b . 14.2b .]) for the streaming sine wave profile so the order of convergence is 
not easy to establish analytically. Also, the timestep is so small that the temporal 
error is dominated by the spatial truncation error. 

4. Test problems. In this section wc compare the different methods from sec- 
tion ([3]) , as they are applied to solve the regularized equation (Eq. [3.1] ) . We consider 
two periodic, one-dimensional test problems: an initial sine wave given by Eq. (|2.ip 
and an initial square wave. We use e = 0.1 unless mentioned otherwise. 

4.1. Smooth initial profile. Figure (14. ip shows profiles at <: = 0.02, beginning 
with the initial profile in Eq. p.l|) . for explicit methods using At = eAa:^/3 Fig ()4.1b : 
consistent with the timestep limit in Eq. [33]) and At = 0.002 Ax (Fig [iHi]'). While 
profiles converge to the physically consistent solution with increasing resolution when 
At cx Ax^, increasing resolution with At cx Ax results in unphysical results as the 
resolution is increased. The convergence of At cx Ax^ and the non-convergence of 
At cx Ax is seen clearly from the plot in Figure (|4.3b .). While Li error decreases for 
At cx Ax^, it increases for number of grid points larger than 16 for At — 0.002Ax. 
Non-convergence starts once the timestep is longer than the limit in Eq. p. 31) , which 
happens for n > 16 for At — 0.002Ax. Although Eq. p.3p is the stability time limit 
on the diffusion term in Eq. p.2p . results do not blow up because diffusion cx 1/e is 
only limited to a small length-scale cx e. 

Quite surprisingly, the numerical experiments show that the implicit methods 
also require a small enough At for convergence. This is probably because of the need 
to accurately resolve diffusion over small length scales (cx e). With e = 0.1, only for 
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Fig. 4.2. The initial profile (Eg. fjlf ) and profiles at t = 0.02 using At = 0.005Ax for 
(a) implicit schemes with e = 0.1 , and (b) the linearized implicit method for e = 0.01 and with 
At = O.OO5A2:, 5 X 10~*Aa;, 5 x 10~^Aa;. Results with the linearized implicit method and the 
nonlinear implicit method look similar. Even for implicit methods, converged results are obtained 
only when At is sufficiently small. 



Ai < 0.005Ax do the implicit methods give non-oscillatory convergent results (we 
have verified this for the number of grid points up to 71 = 262144 — 2^®). Figure 
()4.2b ) shows profiles aX t = 0.02 with implicit methods for the sine wave problem for 
e = 0.1 and At — O.OOSAx. One gets converged results for Ai oc Ax, unlike with the 
explicit methods (it is unclear if it holds as Ax — >■ or if non-convergence appears for 
number of grid points more than the maximum we have tried) . 

Figure (|4.2b ) shows profiles for e — 0.01 but with three different choices for 
Ai, using the linearized implicit method (the nonlinear implicit method gives similar 
results). Using the same Ai that gave converged results for e = 0.1 (At = 0.005Ax), 
one does not get converged results for 71 > 16 (see Fig. |4.3b ]'). Even using a ten times 
smaller timestep (At = 5 x 10^"* Ax; naively one would expect the stable timestep limit 
to scale with e because at extrema, where oscillations originate, tanh[/'/e] w /'/e) 
results in non-convergence for n > 64 (see Fig. 14.3^ ). However, with At = 5 x 
10~^Ax one can see convergence until n = 4096 as seen in Figures (|4.2b ) & (|4.3b ): 
we have verified convergence until n — 32768 = 2^^ for this case. Although the plot 
in Figure (|4.2b ) was done with the linearized implicit method, a similar conclusion 
is reached for the nonlinear implicit method; namely, the nonlinear iterative method 
converges only if At is similar to the value of At which gave converged results with 
the linearized implicit method . 

4.1.1. Numerical Convergence. Figures (|4.3p & (|4.4p show the error (for the 
solution at t = 0.2) for different methods with the initial sine wave test problem (Eq. 
|2.1| ): Figure (|4.3t i) shows errors with explicit methods (using e = 0.1) for At cx Ax^ 
and At c>c Ax; Figure ()4.3b ) shows errors with implicit methods (using e = 0.01) for 
different At scalings; Figure (14. 4p shows errors using e = 0.1, with a stable timestep, 
for all the different methods discussed in this paper. Since the analytic solution 
is unknown, we show the Richardson errors. The Li Richardson error is given by 
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(a) explicit methods with e = 0.1 (b) linearized implicit method with e = 0.01 



Fig. 4.3. Li Richardson error (at t = 0.02 for the initial profile in Eq. }2.1f } as a function 
of number of grid points for (a) explicit methods with e = 0.1 (corresponding to Fig. \4-lV "■'"■'i 
(b) the linearized implicit with e = 0.01 (corresponding to Fig. \4.2\ >]) method using different At. 
Timesteps are chosen to be At = eAx^/S and At = 0.002Ax for the explicit methods, and At ex Ax 
for the implicit schemes. As seen in Figs. ^TTp and converging results are obtained only for 

a sufficiently small At. First- and second-order convergence is indicated by solid lines. 



^"=il/j ~ /il/'^j where fi is the numerical solution for /, with n grid points, at the 
i^^ grid point, and is the interpolation of the solution with 2n grid points at the 
same location. Figure (j4.3|) clearly shows that only for a sufficiently small At do the 
different (both explicit and implicit) methods converge. 

Figure (|4.4p shows that the convergence properties of different methods, using a 
stable timestep, are quite similar; in fact, the errors lie almost on top of each other for 
upwind, linearized implicit, and fully nonlinear methods. All methods show a close 
to second-order convergence in the Li norm. Somewhat surprisingly the Li error 
at the highest resolution is maximum for the semi-implicit method which, unlike 
other schemes, is formally second order in time (however, L2 errors for the semi- 
implicit method were smaller at the maximum resolution). As noted before, since /' 
is discontinuous for the solution, a simple order of convergence does not apply. 

4.1.2. Nature of the converged solution. Now that we have established that 
quite small timesteps are required to obtain converged results, even for implicit meth- 
ods, we will examine the nature of the solutions in more detail. Figure (|4.5p shows the 
time evolution of various quantities (/, v = — tanh(/'/e), vf, and /') for e — 0.0005 
using the linearized implicit method. Figure (|4.5b ) shows the converged profile of 
/ at different times for the sine wave test problem (Eq. 12.11) . using the linearized 
implicit method. The extrema of the sine wave are fiattened quickly and the /' « 
region spreads out in time. The f ^ portion is advected down the gradient at the 
streaming velocity, as expected. The initial sine wave is fully diffused out in a short 
time {t « 0.08). The derivative of the solution changes discontinuously from the flat 
to steep portions (see Fig. 14.5^ . 

Fig. (|4.6p shows the profiles of various quantities (/, v = — tanh(/'/e), vf, and 
/') at t = 0.02 for different es using the linearized implicit method. Timesteps are 
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Fig. 4.4. Li Richardson error as a function of number of grid points at t = 0.02 with different 
methods for the initial profile in Eq. f2. 1}) . The errors for Lax-Wendroff, linearized implicit, and 
nonlinear implicit methods lie almost on top of each other. The solid black lines show scalings as 

and . 

chosen to be short enough that the sohitions are converged for every e. The number 
of grid point is fixed to be 2048. The profiles of different quantities converge as 
e — >■ 0, as required if the regularized equation (Eq. 13. ip is well-posed. Although /' 
changes discontinuously for the converged profile (in regions where solution transitions 
from being flat to steep), somewhat surprisingly v — — tanh(/'/e) is continuous. As 
expected, v = 1 and is down the gradient, in the steep portions of the solution. In 
flat portions, v smoothly changes from -1 to 1, crossing where /' = 0. Although /' 
appears to be zero throughout the flat portion of the solution, it is strictly zero only 
at two points corresponding to the extrema. In flat portion of the solution /' is small 
(but not 0!) and scales with e. The flux of / due to streaming, vf, is continuous, and 
is zero only at points where /' = 0. 

Numerical solution of the regularized equation (Eq. 13. ip is analogous to the 
numerical solution of the Euler/Burger's equations with explicit viscosity. For these 
equations the discontinuous profile at the shock is resolved by viscous spreading over 
a finite scale due to explicit viscosity. Thus, the profiles for various fluid quantities are 
continuous. Similarly, because of regularized diffusion applied where /' ss 0, various 
key quantities (/, v, vf shown in Figs. 14.51 and 14. 6p become continuous. Similar 
smoothing of the regularized solution also occurs for an initially discontinuous profile, 
e.g., the square pulse initial condition discussed in the next section. Fig. ()4.4p shows 
that the solution converges as the resolution is increased for a fixed e. Similarly, 
Fig. (|4.6I) shows convergence for a fixed number of grid points as e — J' 0. Thus, the 
solutions of the regularized equation (Eq. 13. ip converge to a unique solution. 

4.2. Initial square pulse. In addition to the smooth initial profile, we test 
different methods with an initial square pulse. Discontinuous cosmic ray pressure 
profiles can arise in many circumstances, e.g., supernova shocks. The initial profile is 
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Fig. 4.5. Various quantities starting from t = till t = 0.08 (separated by At = O.Olj as a 
function of x for the sine wave test problem (Eq. \2. l\l using e = 0.0005. The linearized implicit 
method with the corresponding stable timestep At = 5 X 10~^ Ax is used. The number of grid points 
is 2048. While f, v = — tanh(/'/e), vf are continuous in space, /' is discontinuous where the 
solution transitions from a fiat to a steep profile. 



given by 

(4.1) /(^,0) = {2 



for 0.4 < x < 0.6, 
otherwise. 



Figure (|4.7p shows the solution of Eq. (|3.ip using exphcit (with A< ~ eAa::^/4, 
consistent with Eq. |3.3| : e = 0.1) and implicit (with At — O.OOSAa;) methods dis- 
cussed in section ([3]) for the square pulse problem (Eq. |4.1| ) . All methods except the 
linearized implicit method give converging results with increasing resolution. Unlike 
with the smooth initial profile, the linearized method fails because it involves taking 
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Fig. 4.6. Various quantities (att = 0.02 ) as a function of x for the sine wave test problem (Eq. 
\2. jl ) using e = 0.1 (blue dashed line), 0.01 (red dotted line), 0.001 (green short-dashed line), & 0.0005 
(solid line). The profiles with different es lie almost on top of each other. The linearized implicit 
method with the corresponding stable timestep (At = 5 X 10~^, 5 X 10~^, 5 X 10~^, 5 X Ax for 
t= 0.1, 0.01, 0.001, 5 X 10~^ , respectively; identical results are obtained for converged solutions 
with other methods) is used for different es. The number of grid points is 2048. Converged profiles 
are obtained as e 0. 



the derivative of the discontinuous square wave; the hnear expansion about /' used 
in Eq. p.7p is not vahd when / is discontinuous. This spurious behavior for higher 
resolutions is eliminated if explicit diffusion (with diffusion coefficient = Aa;; recall 
that the maximum streaming velocity equals unity) is applied. Thus, the linearized 
implicit method (or a semi-implicit version of it) which can be solved exactly using 
a fast tridiagonal solver, is a competitive approach even with discontinuous profiles if 
small explicit diffusion (oc Ax) is applied. 
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Fig. 4.7. The initial profile (Eq. \4-lV o,nd profiles at t = 0.04 for (a) explicit and (h) implicit 
methods at different resolutions with e = 0.1. The linearized implicit method gives unphysical results 
for higher resolutions. However, converged results are obtained if explicit diffusion (with the diffusion 
coefficient = Ax ) is applied in addition, for the linearized implicit method. 



5. Conclusions. We show that the equation governing streaming down the gra- 
dient, although appearing deceptively similar to the advection equation (Eq. I1.2p . 
has an entirely different character. It behaves like a diffusion equation at the extrema 
where /' « (see Eq. I3.2p . and like an advection equation away from extrema. Flux 
moves from the initial maximum to the initial minimum in the sine wave test problem 
(initial maximum is reduced and initial minimum is increased) even after a short time 
(see Fig. |4.5b .]) , while the almost flat maximum and minimum are connected by a 
simple advecting profile. This transport of flux to large distances is a clear sign of the 
non- hyperbolic nature of Eq. (|1.2[) . Similar signs of large scale transport are seen in 
the long tails for the square wave problem in Figs. (|4.7p and (|5.ip . 

The timestep limit for explicit methods on Eq. (jl.2p is severe {At cx Ax^ ; see Eq. 
|2.2| ). However, the regularized equation (with a fixed e) is diffusive and the timestep 
limit for explicit methods scales as At cx Ax^ (Eq. |3.3j ) . The most practical are the 
implicit methods, the timestep limit for which scale as Ax. Although the linearized 
implicit method gives grid scale oscillations for discontinuous profiles, a small explicit 
diffusion cures this problem (see Fig. |4.7b ] V A fully implicit method using GMRES 
does not suffer from such problems and does not require explicit diffusion. Somewhat 
unexpectedly, the timestep limit for getting converged solutions does not scale linearly 
with e; a smaller timestep is required for converging results because of an extremely 
large diffusion (scaling with 1/e) at extrema. It is non-trivial to obtain the timestep 
limit for the convergence of implicit methods (and to establish whether At cx Aa; 
results in a stable scheme as Ax — 0). Traditional linear stability analysis methods 
(e.g., von Neumann stability analysis) are not useful because the streaming equation 
is a highly nonlinear equation even in its most basic form (Eq. |1.2] ). 

Total number of operations with n, the number of grid points in each direction, 
scale as: n'^^'^/e for explicit methods (where d is number of space dimensions); and 
as v}'^'^ / e for implicit methods, assuming that each implicit update scales as n'^ (only 
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Fig. 5.1. Profiles at t = 0.04 for the periodic square wave test problem (Eq. ]4- Ij i evolved with 
the Lax-Wendroff method with At = eAx'^/4) for different box sizes (L), number of grid points (n; 
such that the number of grid points per unit length is the same for all cases) and es:(a) the profiles 
with L = 1, n = 512 are similar for e < 0.1; (b) profiles for L = 16 and n = 8192 (2^^); (c) profiles 
for L = 32 and n = 16384 (2^'^); (d) logarithmic plot for f — 1 as a function of x shows that much 
smaller e is required for convergence with a larger box-size. 



if the matrix is well behaved will Krylov methods converge as n'^; however the exact 
tridiagonal solver will definitely scale as n'^). 

Although the analysis of a simple equation like (|1.2p is interesting in its own right, 
the fluid modeling of cosmic rays involves streaming of cosmic rays down their pressure 
gradient, and is analogous to the simple streaming equation. The generalization to 
multi-dimensions is straightforward, either in directionally split or unsplit fashion. 
Unlike with splitting, unsplit methods do not result in a tridiagonal matrix but can be 
solved by fast Krylov subspace methods. One will need to be careful about calculating 
Vs in the numerical implementation of Eq. (|l.ip in multi-dimensions; interpolation 
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of {B ■ V)pc using slope limiters (as with anisotropic conduction; e.g., [TT]) may be 
required to preserve monotonicity. 

The numerical techniques that we have proposed for solving the streaming equa- 
tion (Eq. fl.2| ) are applicable for different initial (even discontinuous initial profiles 
can be handled; e.g., Fig. |4.7| ') and boundary conditions. Figure (|5.1I) shows that 
e required to obtain converged results scales inversely with the box size, as expected 
(since f'/e ^ f/^L must be ^ 1 for tanh to be a good approximation of the sgn 
function). While e — 0.1 is fine to get profiles similar to those obtained by smaller 
es for L = 1, e < 0.001 is required to obtain the e — ^ solution for L = 16 and 
32. Notice that the e — J' solution consists of a square pulse and extended flat tails. 
Another important point is that the width of the square pulse and the value of / 
in the square pulse at t — 0.04 is the same for different box sizes and different es. 
This implies that the flux leaving the square pulse because of advection and diffusion 
terms in Eq. (|3.2[) is independent of e, box size, and the boundary condition. If the 
pulse was just spreading out at a constant velocity (instead of having the extended 
tails seen in Fig. |5.1| ). the value of / at t = 0.04 within the square pulse would be 
/ = 1 + 0.2/(0.2 + 2 X 0.04) = 1.714, much larger than the value of 1.425 seen in the 
converged solutions with streaming. The average rate of decrease for / within the 
square pulse is almost double for the streaming case with tails as compared to a naive 
advection solution. Similar large scale transport of flux is seen in the sine wave test 
problem in Figure (|4.5b ). The relative proportion of fluxes in the advecting square 
wave and the long tails depends on the height of the square pulse; higher pulse height 
leads to a larger flux in the advecting pulse as compared to the extended tails. Thus 
realizing the diffusive nature of Eq. (|1.2|) can be astrophysically important, e.g., to 
accurately obtain the cosmic ray loss rate from the Galaxy and leakage of cosmic rays 
upstream of supernova shocks. 

Comparison of the fluid model (jl.ll) with particle simulations of cosmic rays is 
required to test if the cosmic ray streaming model applies for a particular situation 
(e.g., the tails clearly cannot extend to space scales larger than the light crossing time!) 
but once the validity of the streaming model for the cosmic ray fluid is established 
(probably true for Galactic scales), our inferences are applicable. In addition to the 
implicit diffusive nature of the streaming equation, diffusion of cosmic rays because of 
scattering by Alfven waves (ignored in this paper) can result in additional spreading of 
cosmic rays; the cosmic ray diffusion time is typically much longer than the streaming 
time for typical Galaxy and galaxy-cluster parameters (e.g., [7]). Streaming of cosmic 
rays at the local Alfven speed is extremely important for cosmic ray conflnement 
in the Galaxy for millions of years. Cosmic rays can heat the interstellar and the 
intracluster medium (due to the \vs ■ Vpcrl term on the right hand side of Eq. |1.1| or 
because of convection driven by them) where cosmic ray pressure (compared to the 
thermal plasma pressure) is signiflcant (e.g., [3l[12]). Thus, the methods developed in 
this paper will be useful in modeling cosmic rays interacting with the thermal plasma 
in various astrophysical contexts. 

Solving the regularized form (Eq. 13.11) of the streaming equation (Eq. 11.21) results 
in a unique, converged solution, provided that the timestep is small enough. Regular- 
ization is analogous to introducing explicit viscosity in solving Euler/Burger's equa- 
tions and resolving the viscous length scales at the shocks. With regularization the 
solutions are approximated very accurately, the advection velocity (v = — tanh(/'/e)) 
and streaming fluxes (vf) are continuous, and converged results are obtained in the 
limit e -> and Aa; ^ (e.g., see Figs. 14.41 fc lT^ . 
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More work is clearly needed to understand the nature (both numerically and an- 
alytically) of the streaming equation (Eq. II. 2p . For example, like the Euler equations 
which can result in discontinuous profiles starting with a smooth initial condition, a 
smooth initial profile for the streaming equation can give rise to discontinuous first 
derivatives (e.g., /' changes discontinuously from w within the flattened sine wave 
[e.g.. Fig. 14. 5| to a constant value outside of it). There is significant transport over 
large scales (e.g., tails in Fig. |5.1| ) because of the parabolic nature of the streaming 
equation at extrema. 
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